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ABSTRACT 

An implicit Lagrangian hydrodynamics code for general relativistic spherical 
collapse is presented. This scheme is based on an approximate linearized 
Riemann solver (Roe type scheme) and needs no artificial viscosity. This 
code is aimed especially at the calculation of the late phase of collapse-driven 
supernovae and the nascent neutron star, where there is a remarkable contrast 
between the dynamical time scale of the proto-neutron star and the diffusion 
time scale of neutrinos, without such severe limitation of the Courant condition 
at the center of the neutron star. Several standard test calculations have been 
done and their results show (1) this code captures the shock wave accurately, 
though some erroneous jumps of specific internal energy are found at the contact 
discontinuity in the shock tube problems. (2) The scheme shows no instability 
even if we choose the Courant number larger than 1. (3) However, the Courant 
number should be kept below ~ 0.2 at the shock position so that the shock 
can be resolved with a few meshes. (4) The scheme reproduces the well known 
analytic solutions to the point blast explosion, the gravitational collapse of 
the uniform gas with 7 = 4/3 and the general relativistic collapse of uniform 
dust. Two other adiabatic simulations have also been done in order to test the 
performance of the code in the context of the collapse-driven supernovae. It is 
found that the time step can be extended far beyond the Courant limitation at 
the center of the neutron star. The details of the scheme and the results of these 
test calculations are discussed. 

Subject headings: Supernovae — general relativity — hydrodynamics - 
numerical simulation 
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1. INTRODUCTION 

Since the pioneering work by Colgate and White (1966), detailed studies of the 
dynamics of collapse-driven supernovae have been done mainly with numerical simulations 
( Bcthc 1990| , Miillcr 1991 , and references therein). One of the difficulties in so doing is 



the remarkable contrast between the dynamical time scale of the nascent neutron star 
(<z 1msec) and the diffusion time scale of neutrinos (~ lsec). If we want to simulate 
the whole scenario of the collapse-driven supernovae, these time scales should be treated 
simultaneously. Since the typical time scale of weak interactions (~ 10~ 7 sec) is much 
smaller than the dynamical time scale, the neutrino transfer has been treated in implicit 
ways in general. The hydrodynamics, however, has been calculated chiefly by explicit 
schemes, so that the time steps are restricted by the Courant condition and a large number 
of integration steps are required to simulate the late stage of collapse-driven supernovae. 
Two approaches have been conceived to overcome such severely limited time steps. The 
first one is to implement an algorithm of individual time steps in the explicit difference 
scheme ( Powers fc Wilson 1991| ), in which the time steps at different positions are allowed 



to be different from one another and are determined from the local Courant condition. The 
advantage of this approach is that it is easily applied to multi-dimensional simulations. On 
the other hand, the implicit difference scheme is another possibility, since the time step is 
not restricted by the Courant condition, though the required accuracy of calculation will 
limit the time step. The disadvantage of this approach is that the number of operations per 
each time step becomes much larger than in the corresponding explicit scheme, particularly 
in multi-dimensional simulations. Hence, it is of no use unless we can take time steps large 
enough to make up for its larger operation number. However, note that the nascent neutron 
star, whose dynamical time scale limits the time step in the explicit schemes, evolves in a 
nearly hydrostatic manner so that it might be possible to take large enough time steps in 
this case. Note also that an advantage of implicit hydrodynamic schemes in the simulation 
of collapse-drive supernovae is that the neutrino transfer schemes are generally coded in 
implicit fashions so that the hydrodynamic parts are easily incorporated into the neutrino 
parts. 

The purpose of this paper is to provide such an implicit scheme from the latter stand 
point. Such attempts have already been made by Schinder et al. (1988) and Swesty (1995). 
Their schemes use almost the same equation set as that of May & White (1967), which 
has been frequently used in explicit general relativistic calculations (note, however, that 
Schinder et al. extended it to polar-slicing coordinates). Their codes also need artificial 
viscosities. On the other hand, for the present scheme we tried to apply a so-called 
Godunov-type differencing method ( [Hirsch 1990| and references therein), which has been 



extensively utilized for the multi-dimensional hydrodynamic simulations these days. Since 
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we have evaluated the pressure and velocity at a cell interface using the solution of the 
linearized advection equations, this code is an extension of Roe's scheme. Thanks to the 
numerical diffusion induced implicitly by finite differencing, this code needs no artificial 
viscosity like most of the recent elaborate Euler schemes. In addition, it is quite robust. 
In fact, as shown later, the calculation does not collapse even if we take time steps much 
larger than the Courant limit, (though the results become less accurate). 

This paper is not aiming at the general application of Roe's scheme to general 
relativistic hydrodynamics QEulderink fc Mellema 1995| , [Romero et al. 1995|) . The 
coordinates and the metric we use here are fully reduced ones under the assumption of 
spherical symmetry so that the degree of time slicing does not remain at all. The equation 
set is also so chosen that they look similar to the nonrelativistic counter parts in these 
coordinates. Instead, the final goal of this project is to study the neutrino physics in 
collapse-driven supernovae by combining this code with sophisticated neutrino transfer 
codes. Recently, many researchers in this field have considered that the key factors of 
successful explosions are neutrino and multi-dimensional hydrodynamical effects such as 
convection ( [Her ant et al. 1994 , Burrows et al. 1995 , Janka & Mullcr 1995[ [Sliimizu et al 



|1994|) . It is true that the spherical simulations cannot treat such multi-dimensional effects 
properly, but we think that they can play a complementary role to mult i- dimensional 
simulations since mult i- dimensional simulations use a more approximate treatment of 
neutrino transport than in spherical calculations. At present we are trying to incorporate 
this hydrodynamical code into the multi-group flux- limited diffusion code ( [Suzuki 1990|) . 



The coding of the Boltzmann solver has also been undertaken. This paper is the first step 
of the project. 

The organization of this paper is as follows. The details of the formulations and the 
schemes are described in the next section. The results of the test calculations are shown in 
the section 3. Concluding remarks are provided in the final section. 



1.1. Basic Equations 

Spherical symmetry of the system is assumed and the following form of the metric is 
used: 

In the above formula c and G are the velocity of light and the gravitational constant, 
respectively, which are taken to be unity in the following equations, t is the coordinate 
time and m is the baryon mass coordinate which is related to the circumference radius r 
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through the conservation law of baryon mass as described below. This form of the metric is 
so simple that it has been frequently adopted in the spherically symmetric simulations thus 
far. However, this way of time slicing is not suitable for the study of black hole formation 
(see |Schinder et al. 1988| ). Since our purpose is to investigate the supernova dynamics 
leading to proto-neutron star formation, that is not a serious problem. The basic equations 
consist of the Einstein equations obtained from the above metric: 

G» v = 8ttT^ (2) 

and the Euler equations: 

V U T^ = (3) 

= { Pb {l + e) + p}uV - pg» v (4) 

with the baryon number conservation equation: 

V> 6 ^) = (5) 

and the evolution equation of electron fraction (see below). Here G^ u and T^ v are the 
Einstein tensor and the energy-momentum tensor, respectively. p\> is the baryon mass 
density and m m is the four velocity of the matter, e and p are the specific internal energy 
density and the matter pressure, respectively. g^ v is the inverse of the metric tensor g^. In 
addition to these equations, the neutrino transfer equations should be implemented finally, 
but they are dropped in this paper. For this reason, the source term which evolves the 
electron fraction is not included in the present calculations. First I will write down the 
hydrodynamic equations with those representing the baryon number conservation and the 
evolution of electron fraction. They are formulated as follows: 

e'^^r = -t^ 4vrr 2 m — - 6 

dt Tdm K ' Tr K J 

e -*?E = _E 4 7rr 2 + ^) - ™ 

dt h \ dm Anr 2 J r 



-T 47rr I + irhi - — -4irr(p+p v ) (7) 



-<b de 1 f d i a 2m ^r 2 rF u \ 

'dt = -f{ P ^ rU) ~ P ^—)- rQ 

-J^ 2 ^— (p + p v ) - -rUq + y -^rrF v - tQ (9) 
,dY e r(5f v \ d 3 p 

1* = ^ T J\^UT (10) 



<9(|7rr 3 
dm 



Tr . (11) 
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In the above equations, r = — is the inverse of the baryon mass density, and 

Pb 

Ot 

U = e - * — = D t r is the radial fluid velocity. T is the general relativistic gamma factor 
at 

which is defined below (equation fllED) and is solved simultaneously, h is the specific 
enthalpy which is also defined below (equation (|TPD) and is treated as the dependent 
variable to be solved, m (r) is the gravitational mass inside the radius of r and should be 
distinguished from the baryon mass coordinate m. The second representation of the energy 
conservation @ looks peculiar. This form of the equation declares explicitly that we get the 
internal energy density by solving the equation of the total energy conservation of the fluid 
(in the Newtonian limit) and then extracting the kinetic energy and the gravitational energy 
from it rather than using the simpler form (§) which is just the first law of thermodynamics. 
Both forms of the energy equation are tried in the shock tube calculations and it is found 
that the latter form (^) better reproduces the Rankine-Hugoniot relation. This is a well 
known fact in the Eulerian numerical simulations. Note that the artificial viscosity is not 
necessary in both cases for the finite difference scheme described below. F u , p u , q, Q and /„ 
are quantities related with neutrinos. They are defined as follows: 

Tf = p b E u u a uf > + F?u f> + u a FP + pf (12) 
Q a ee T^ ;Q = (e-*Q,e- A g,0,0) (13) 




P^ ■ (14) 
coll. P 

Tf is the energy- momentum tensor of neutrinos and pb is the baryon mass density. E v , 
F° and p°f are, respectively, the energy density, the flux vector and the stress tensor of 
neutrinos. The right hand sides of equations fllQ|) and ([H]) are the phase space integrals of 
the collision term in the neutrino Boltzmann equations. Although these terms are included 
in the above equations for completeness, they are omitted in the present calculations. In 
particular, the right hand side of equation (pi]) is set to be zero, so that Y e = const, in all 
of the computations below. Equation (|TTD relates the radius and the baryon mass inside it 
and is obtained from the combination of the baryon number conservation equation (|5[) and 
the definition of F below (equation (0)). This equation is used to determine r(m) of the 
next step instead of D t r = U, which has been often used so far. It should be noted that the 
above equations are very similar to the Newtonian counterparts. As a result, the difference 
scheme developed for the Newtonian system can be applied to the above general relativistic 
system. The rest of the basic equations specify the metric components X(t, m) and <p(t,m), 
the gamma factor T(t, m), the gravitational mass m (t, m) and the specific enthalpy h(t, m). 
These equations are: 

\ 1 or , . 

e = (15) 
Tdm K J 
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T 2 = 1 + u 2 16 

r 

h d(j) = _ r dp_ _ r 2 g 

dm dm 47rr 2 

f) m f) A 

°— = -(p + p^i-nr^-Arrr'eTF^ (18) 

h = l + e + pr . (19) 

As is clear, these are not evolutionary equations except for equation (|18|), but are constraint 
equations which must be satisfied at each time step. Since the present scheme is written in 
the implicit form, these constraints are always satisfied within the error for convergence of 
iteration. The equation for the gravitational mass m can be written in the constraint form 
as well: 

However, equation (|18D is adopted here rather than (|20|), since it explicitly guarantees the 
total energy conservation m{t, R) = const., where R is the stellar core radius, as long as the 
energy loss due to the emission of neutrinos is disregarded. The constraint equation ( pUD is 
used for the check of the accuracy of calculations. The terms concerning the interactions 
of neutrinos are included in the above equations also, but are dropped in the present 
simulations again. Equation (^) is nothing but the definition of specific enthalpy. Usually, 
this expression is expanded in the preceding equations. However, it is treated as one of the 
dependent variables to be solved simultaneously in this code, since it avoids the repeated 
appearance of the derivatives of h in the matrix elements. This treatment means that the 
specific enthalpy does not satisfy the above definition strictly during the iteration process. 
This did not cause any serious problem in the test calculations shown below. Moreover it 
increases CPU time very little. It should be noticed that the most time consuming part in 
the radiative hydrodynamics code is the neutrino transfer subroutine. In order to integrate 
equation (|I7|) we must specify the boundary condition for 0. It is done so that the metric 
in the stellar core coincides at the surface of the stellar core with the Schwarzschild metric 
for a spherically symmetric vacuum: 




2m„ 



(21) 



where the subscript "s" means the variables are evaluated at the surface of the stellar core. 
Strictly speaking, the metric outside the stellar core is not the Schwarzschild metric because 
neutrinos exist there. For simplicity, however, the above approximation is adopted in the 
following calculations. To summarize, the basic equations to be solved are equations (pf) 
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through (|ll|) and (15) through (0). The latter form (|9|) is used for the energy conservation 
equation as stated above. The boundary condition is given as equation (^1|). These 
equations are differenced in the implicit method described in the next section. The resulting 
equations are the simultaneous equations consisting of 10 x N (mesh number) equations. 
The matrix to be inversed has a block-tridiagonal form. 



1.2. The Implicit Differencing Scheme 

All the dependent variables except r, <fi and m are defined at the mesh centers as 
in most of the recent Eulerian schemes. The values at the cell interfaces for pressure 
and velocity are calculated by using the solutions of the linearized Riemann problems as 
described below, r, and m are placed at the edges of meshes in order to simplify the 
integrals (equation (fPTD) and the implementation of the boundary condition (equation fl2T|)). 
The superscripts and the subscripts attached to the variables in the following difference 
equations are representing the time step and the mesh number, respectively. It is assumed 
that the mesh centers are specified by the lower case subscript like "i", on the other hand 
to the interfaces attached are the uppercase subscript like "I", and that all these subscripts 
are integers. It follows that i-th mesh center {i = 1 ~ N) is located between the (/ — l)-th 
and J-th interfaces (/ = 1 ~ N). However, there is one exception. 0/s are defined on 
the (I — l)-th interfaces. This is because the boundary condition of </> (equation (pT|) ) is 
specified on the surface of the iron core (iV-th interface). This is not the case for r and m, 
for which the inner boundary conditions: r c = and m c = are assumed implicitly and we 
don't have to solve them. 

Paying attention to these assumptions, we will first see the finite differenced equations 
regarding the z-th mesh and 7-th interface (but (/ — l)-th interface only for 0): 



= Am * C -<£»I' 
4ttA£ 




+ A ~ ^ZnU + Am,<r>r* < P >i 2 (23) 



47r rfr- 
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K +1 " (l + + rj' 



.n+l n+l 
Pi 



(24) 

(25) 
(26) 
(27) 

(28) 

(29) 
(30) 
(31) 



In the above equations it is supposed that the time step "n" corresponds to the present time, 
where all the values of the dependent variables are known. The above finite differenced 
equations form nonlinearly coupled simultaneous equations. The Newton-Raphson scheme 
is used to solve these equations. Some averaging procedures are abbreviated like < • • • >, 
<C • • • and 7TT . Their definitions are as follows: 



X n+1 + X? 



i V n+1 _l_ V n +^ J_ Vn _l_ V n 



:zrf + (zf_ 1 y 



(32) 
(33) 

(34) 
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where X is assumed to be defined at the cell center or at the cell interface while Y is 
supposed to be defined at the cell edge. Z is actually the radius r. Note that we distinguish 
two velocities U™ +1 and and two pressures and p™ +1 , the latter of which are 
determined from the linearized Riemann problems which are evaluated at the time step 
"n + 1", where the values of the dependent variables should be solved. That is, we consider 
the Riemann problems at each cell interface "I" at each time step "n + 1", the left and 
right constant states of which are given by (r™ +1 , and (r^ 1 , U^ 1 , p™+i), 

respectively. We solve these problems approximately by linearizing the advection equations 
as described below. It should be noted again that these constant states are evaluated 
at the time step "n + 1", so that we cannot know them a priori. Hence we must solve 
these Riemann problems at each iteration step of Newton-Raphson methods. The physical 
meaning of the interface values obtained in this way is rather obscure compared with the 
counterparts in explicit schemes, since the exact solution of the Riemann problem describes 
the interaction of the nonlinear wave after the time of the initial condition (the time step 
"n + 1"). The interpretation here is that they approximate the states after the interactions 
of nonlinear waves during the time steps "n" to "n + 1". Practically, these values are given 
as follows. If we pay attention only to the advection terms, then we get 




V 



— 4nr 2 

r 

o 



P e A 2 

—7 Airr 

t r 



where 7 is the adiabatic index and is defined as: 





e^r 



\ 



h 




(35) 



7 



d\\ip s 
din pi,, 



[s is the entropy per baryon) 



(36) 



In equation (0), the dependent variable is changed from e to p, for two reasons. Firstly, 
the eigen vectors are very simple in this representation: 
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(37) 



(38) 



2Tc x 



(39) 



V 2 1P J 



where A's are the eigen values, and r's and Z's are the corresponding right and left eigen 
vectors. In this equations, c s is the sound velocity defined as: 



V 
p b h 



(40) 



As in the Newtonian hydrodynamics, there are three modes, two of which correspond to the 
left- and right-going sound waves and the rest of which represents the contact discontinuity. 
In this representation, it is clear that pressure and velocity are continuous across the contact 
discontinuity. This feature is taken over to the solution of the linearized Riemann problems 
when these eigenvectors are used to solve them. It is not the case when other dependent 
variables are chosen. This is the second reason to use these variables. In order to linearize 
equation ( p5|) we just replace the Jacobian matrix by some constant matrix at each cell 
interface. In this paper, arithmetically averaged p bm = (p b L + Pwj)/2, U m = (Ui + Ur)/2 
and e m = (el + £r)/2 are used to evaluate the constant matrix, where the subscripts L 
and R mean that they are evaluated at the left and right constant states of each Riemann 
problem, respectively. The eigen values and the corresponding left and right eigen vectors 
given above are also evaluated in terms of these averaged variables at each cell interface. 
The linearized advection equations thus obtained can be solved analytically, ui and pi 
evaluated in this way prove to be as follows: 



Uj 



Pi 



-{U R + U L ) 



^ImP; 
^mPm 



^(Pr + Pl)- 2c 



-{Pr-Pl) 

n 

{u R ~ U L ) 



(41) 
(42) 



where the subscript m means they should be evaluated by using the above averaged 
variables. It should be noted that p m is different from the arithmetic mean of pl and pr. 
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As is clear, both ui and pi are the arithmetic mean of the left- and right-state values with 
the correction terms which are proportional to {p^ — Pl) an d (ur — ul), respectively. These 
correction terms serve as a numerical diffusion so that no artificial viscosity is necessary. 
The above representations of the velocity and the pressure at the cell interface are so simple 
that it is quite easy to evaluate their derivatives with respect to the dependent variables in 
the course of the Newton-Raphson scheme. 

In the present code, in order to achieve spatial second order accuracy, the piecewise 
linear distributions of dependent variables pb, U , e and T are introduced. It is well known 
that some limiting procedure is necessary to maintain the monotonicity of the numerical 
solution. In this paper, we use the simplest slope prescription: if s t and s r denote the 
left- and right-hand difference quotients of the above dependent variables, the slopes S are 
determined by the formula 



for |s|| < \s l r \ and s\ ■ > 

for \s\\ > \sl\ and s\ ■ s l r > (43) 



otherwise 



where the superscript i distinguishes the dependent variables p&, U, e and Y. This 
prescription will need some improvement in the future. 



1.3. Numerical Implementation of the Implicit Scheme 

The above finite differenced equations (|22D through ( j31~l) with the definitions of interface 
values ( |4~TD and (|42|) form the system of the nonlinearly coupled simultaneous equations 
with respect to the dependent variables r, U , e, Y e , r, A, T, 0, m and h. In order to solve 
these equations, the standard Newton-Raphson iteration scheme is utilized for the linearized 
equations. That is, if we define the dependent variables vector as 

X =(•■•, n, Ui, Ei, Y ei , n, Xi, Ti, 7 , rhj, hi, • ■ •)* , (44) 

where the subscripts i and I mean the i— th mesh and interface values, and denote all the 
equations in the vector representations as 

F(X) = , (45) 

then the linearized equations become 

F(X) + ^P--5X = , (46) 
dX 
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where X is some guess for the correct X. This linearized equation fl46f) is solved repeatedly 
while improving the guess by 

Xnew = X + 5X (47) 

until the error \SX\ becomes less than some small value which is set to 10~ 5 in the following 

dF(X) 

calculations. As is clear, the Jacobian matrix — — turns out to be a block-tridiagonal 

dX 

matrix whose block size is 10 x 10. At each step of the above iteration, the linearized 
equations ( pf) are solved by the standard Gauss elimination scheme, though other elaborate 
schemes can be applied to this type of matrix. Note, however, that the matrix for the 
hydrodynamic sector is very small compared with that for neutrinos. The implementation 
of such an efficient solver will be done at the next step, where the neutrino transfer code 
will be combined. 

As stated above, the formulae for the velocity and the pressure at the cell interface 
(equation ( |4"I]) are so simple that we can calculate their derivatives with regard to the 
dependent variables analytically. In so doing, however, the contributions from the sloping 
prescription are dropped, which does not cause serious problems for convergence of the 
Newton-Raphson iteration process. The thermodynamic derivatives are obtained by 
differencing the EOS table numerically, except for the adiabatic index 7, which is stored in 
the table. 



1.4. Time-Step Control 

The time step is determined so that the change of any dependent variables should 
not exceed an appropriate upper limit, which is chosen to be 2% in most of the following 
calculations, at each time step. In order to save the computation time, however, we correct 
the next time step and don't repeat the same time step even though the above criterion is 
not satisfied at some time step. In fact, if the change of some variable exceeds the above 
criterion, the subsequent time step is decreased by a factor of 0.95. On the other hand, if 
none of the variations of the dependent variables exceeds the upper limit, then we multiply 
the next time step by 1.05. The time variation of the Courant number which is defined as: 

Courant number = —. r- (48) 

\Cs+\U\ 



is monitored for all the calculations and is shown for some test calculations below. 
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2. TEST CALCULATIONS 

In the following, some results for the representative test problems are shown in order 
to clarify the performance of the present hydrodynamical code. Among them there are 
shock tube problems both for non-relativistic and special relativistic cases, Sedov's point 
explosion problem in the uniform matter, Oppenheimer-Snyder's dust collapse problem and 
the non-relativistic self-similar collapse of a 7 = 4/3 polytrope. Moreover, the results of the 
hydrostatic calculations and the simple adiabatic collapse simulations are demonstrated in 
order to judge whether the present code can be applied to the simulation for the late time 
stage of the collapse-driven supernovae. In order to compare the results, the parameter sets 
adopted here for the first five problems are almost the same as those studied by Swesty. In 
the course of these simulations no special tuning to the specific problem is done except for 
the boundary conditions, the mesh system and the equation of state. 



2.1. The Non-Relativistic Shock Tube Problem 



The most standard parameter set for this problem is that proposed by Sod (1978): 



PL 
PL 

U L 



1 
1 




PR 
PR 
Ur 



0.125 

0.1 





(49) 



where the subscripts L and R mean the left- and right-state values and CGS units are 
assumed. The equation of state for this problem is p = (7 — l)p&£ with 7 = 4/3 and 5/3. 

Since the present code assumes spherically symmetric geometry, we formulate the 
problem, which is originally for plane symmetric geometry, over a thin shell whose curvature 
is safely ignored. Following Swesty, we take R = 10 4 cm and 5R = 2cm for the radius and 
the thickness of the thin shell, respectively, and we define x = R — r and specify the region 
— 1 < x < as the left state. 100 uniform meshes are used. Although the calculations 
have been done both for 7 = 4/3 and 5/3, only the results for 7 = 5/3 are described 
since the features of the numerical solutions are nearly the same both qualitatively and 
quantitatively. 

In figure 1 the representative distributions of p&, U, p and e for 7 = 5/3 are shown 
with the corresponding analytic solutions. The time is 0.675sec (400 time steps). As can be 
seen, the shock is resolved by 3 meshes in the present case. Moreover, no oscillations occur 
behind the shock and the rarefaction wave. The Rankine-Hugoniot condition is strictly 
satisfied without any artificial viscosity. However, we can see the small jump of the specific 
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internal energy density ahead of the contact surface. This is not an oscillation but an 
erroneous contact surface. This can be understood from the smoothness of the velocity and 
the pressure at the contact surface. The approximate solution of the Riemann problem is 
least accurate at the first time step, when the largest discontinuities exist. The observed 
incorrect jump of the specific internal energy density at the contact surface is the result of 
the inaccurate approximate solution at the first time step. Since this scheme is Lagrangian, 
the contact discontinuity is never smeared out at all once it is formed. However, this is not 
a serious problem in the simulation of the realistic core collapse, since even the contact 
surface between the iron core and the mantle spreads over a considerable distance and is 
not a strict discontinuity. 

The time steps for these calculations are controlled so that the maximum variation 
of any dependent variables does not exceed 2%. As a result, the Courant number defined 
earlier is kept as large as 0.3 through the calculations. The calculation which allows the 
maximum variation of 15% has also been done in order to see the effect of time step control. 
The result shows the tendency to smear out the shock and the rarefaction wave. However, 
the effect is not so large in the non-relativistic case. This issue will be raised again in the 
next section. 



2.2. The Special Relativistic Shock Tube Problem 



Next shown are the numerical solutions for the following Riemann problem: 



Pl = 0.2 

p L = 2 x 10 18 

U L = 



Pr = 0.1 

p R = lx 10 18 

U R = 



(50) 



where the units are again CGS. This parameter set is the same as that used by Swesty. We 
use 100 zones with an equal spatial resolution also in this case. The equation of state is 
again p = (7 - l)p b e. 

Figure 2 shows the distributions of density, velocity, pressure and specific internal 
energy density at the time of 3.03 x 10 _10 sec (490 numerical steps). 7 = 5/3 also in 
this case. As can be seen clearly, the Rankine-Hugoniot condition at the shock wave is 
satisfied very well. No oscillation occurs in this case, either, except for the one mesh just 
ahead of the contact discontinuity. The cause of this erroneous jump is the same as for 
the non-relativistic shock tube. Since the initial discontinuity of the density is smaller 
in this case than in the previous case, the resulting error becomes smaller. One thing to 
be mentioned here is that the upper limit of the variation of the dependent variables is 
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chosen to be 0.5% in this calculation so that the shock can be resolved in a few zones. The 
resultant Courant number is ~ 0.2. In order to see the effect of time steps on the resolution 
of the numerical solution, we have done two other simulations varying the upper limit. In 
figure 3 we show the results for the upper limit of 2%. The Courant number is ~ 2 through 
the integration. As can be clearly seen, both the shock wave and the rarefaction wave are 
smeared out considerably. This trend is much more remarkable in figure 4, which shows 
the results for an upper limit of 10% and the Courant number is as large as 8 at the end 
of the calculation. It should be also noted that in this calculation the Courant number is 
set to be 3 from the beginning. From these results it is concluded that although the present 
implicit scheme is stable for the Courant number > 1, in order to keep the resolution of 
the solution the maximum variation of the dependent variables must be chosen so that 
the Courant number should be ~ 0.2 at the shock wave. However, note that this does not 
mean that this code cannot not be applied to the realistic simulation of the collapse-driven 
supernovae, since what we should overcome is not the Courant limit at the shock wave but 
the limit at the nascent neutron star. This issue will be discussed again in the section of 
adiabatic collapse calculation. 

2.3. Sedov's Point Explosion in a Uniform Medium 

The dimensional parameters which characterize the problem are the density of the 
uniform matter po and the total energy of explosion -E^ot- The problem, however, can be 
formulated in the completely non-dimensional form. Hence the solutions with different 
initial parameters are equivalent to each other. We choose the following parameters for 
simplicity: 

Po = 1 

< E tot = 1 , (51) 

^ 7 =5/3 

where the units are CGS. The computational region is the sphere with the radius of 1, 
and is covered by 100 equal baryon mass zones. The equation of state is p = (7 — l)pft£. 
The calculation is initiated by depositing the total energy in the central mesh. For 
numerical convenience the uniform ambient matter has the specific internal energy density 
of ~ 10~ 3 cm 2 /sec 2 . This is justified since the ambient specific internal energy density is 
negligibly small compared with that of shocked matter. The results are shown in figure 5. 
The shock is resolved by 2 ~ 3 meshes also in this case, and the analytic solution is 
reproduced quite well in most of the computational region. In the central ~ 10 meshes, 
however, the numerical solution is not very good. This is because the mesh for this problem 
is uniform in baryon mass and the spatial resolution is not sufficient for the central region. 
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2.4. The Self-Similar Collapse of a 7 = 4/3 Polytrope 

It is well known that a non-relativistic gas with an adiabatic index of 7 = 4/3 is 
neutrally stable against gravitational collapse. Goldreich & Weber (1980) found a solution 
of a self-similarly collapsing gas with 7 = 4/3. Yahil (1983) generalized this solution to 
the collapse for different 7's. These self-similar solutions approximate the homologous 
collapse of the inner core of the collapse-driven supernovae. Hence it is indispensable to 
the numerical schemes for the collapse-driven supernovae that they can reproduce this 
self-similar collapse. 

The initial condition is made as follows: assuming the polytropic equation of state 
p = KpJ with 7 = 4/3, we first solve the Lane-Emden equation. Here we take K as the 
value for the relativistically degenerate gas with electron fraction Y e = 0.45. The central 
density is set to be 1.0 x 10 8 g/cm 3 . Then we take out the inner IMq and reduce the 
internal energy uniformly so that the total energy becomes 0. This reduction factor is 
nearly equal to the pressure deficit d defined by Goldreich- Weber as 



also in this case. It is noted that the equation of state for simulation is ideal gas like 
(p — (7 — l)pbs) with 7 = 4/3, not adiabatic (p = KpJ). This is because we want to see 
how well the constancy of entropy is maintained in the course of the numerical integration. 

The results are shown in figure 6, where the time variation is presented for density, 
velocity, pressure and the quantity proportional to entropy: p/pl- In this figure, the solid 
lines are just connecting the data for the same time step. When we pay attention to the 
velocity distribution, it is clear that most parts of the polytrope collapse homologously, 
that is, the infalling velocity is proportional to the radius. Although it is not the case 
for the outer few zones, this is due to the boundary condition which is set by hand and 
is not consistent with the actual dynamics. In the actual simulation of the supernovae, 
this problem concerning the outer boundary condition should be avoided by locating the 
boundary at a sufficiently distant position. On the other hand, it is also found that the 
entropy distribution remains almost constant for most parts of the polytrope through the 
evolution. The exception is the outer few zones also in this case. Ignoring these outer parts, 
we can find the maximum error of the entropy is about 8% at the innermost zone and about 
5% for the other zones. The calculation is terminated when the central density reaches 



2 




(52) 



where M-^ c and M are the homologous core mass which is set to be IMq in this model 
and the mass of the polytrope without pressure deficit, respectively. 

The computational region (IMq) is covered with 100 equal baryon mass zones 
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4.85 x 10 12 g/cm 3 after 600 time steps. The relativistic correction is goo ~ 0.97 at that time. 



2.5. Oppenheimer-Snyder's Dust Collapse Problem 

The general relativistic collapse of uniform pressureless matter was analytically solved 
by Oppenheimer & Snyder. When we assume the matter is initially (t = 0) at rest 
with the uniform density pbo and the radius r , then the subsequent motion is described 
parametrically as follows: 

r = -^(l + cos?7) (53) 



* = ^V^^ (r?+Sinr?) ' (54) 

where r\ is varied from to 7r. It should be mentioned that the above t is the proper time for 
each mass element and should be distinguished from the coordinate time t. The boundary 
condition for the coordinate time t is not altered in this simulation, they are related to each 
other by 

ar e - ■ < 55 > 

The time coordinate in figure 7 is the coordinate time t which is obtained by numerically 
integrating equation (|53p. However, the difference is not significant except for the final few 
msec. Note also that in this time coordinate it takes infinite time for the mass shell to reach 
the event horizon. However, it is not a serious problem, as stated earlier, since the main 
purpose of this code is to study the process to produce a neutron star. 

Following Swesty, we take 2Mq cloud with a uniform density of pbo = 10 8 g/cm 3 for the 
initial model. Since the present code cannot handle zero-pressure, we use the polytropic 
equation of state p = Kp\ with 7 = 4/3 and negligibly small K. This initial configuration 
is covered with 160 equal baryon mass meshes. The calculation is terminated after 3000 
time steps when the central density reaches 1.72 x 10 14 g/cm 3 . 

In figure 7, the motions of the representative mass elements are shown. The solid 
lines are exact solutions in this case. As shown in this figure, the coincidence between the 
numerical solution and the exact solution is satisfactory. The slight difference (less than 1% 
in time) appears only for the last few msec which is not clear in figure 7. This is partially 
due to the error of the integration of the coordinate time (equation (|53"D). The above result 
along with the next one clearly shows that the present code can treat the general relativistic 
gravity correctly. 
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2.6. The Hydrostatic Calculation 

This and the next section are devoted to demonstrate the performance of this code in 
the context of the collapse-driven supernovae, that is, it can calculate the hydrodynamics 
far beyond the dynamical time scale after the proto-neutron star is formed in the central 
region. This is the crucial point for the purpose of this paper. 

In this section, we take a hydrostatic configuration for the initial condition for the 
dynamical calculation, and see whether this initial state will be maintained or not. For that 
purpose, we first solve the Tolman-Oppenheimer-Volkov equations: 

dp m(r)+Anr z p 

-j- = -(P + P) 7 — ^TT\ ( 56 ) 



d<j> 

dr 




-2 



2 m (r) 



(57) 



r 

P = Pb{l + e) (58) 

m (r) = r p(r')Airr' 2 dr' , (59) 
Jo 

where the notation is the same as in the dynamical equations. The boundary condition for 
(f> is consistent with the dynamical counterpart: 

f .- = [l_^£) . (60) 



The subscript s in this equation means the core surface values. The equation of state is 



assumed to be polytropic p = KpJ with K = 1.97 x 10 3 in the CGS units and 7 = 2.5. 



The central density is taken to be pb c = 4 x 10 14 g/cm 3 . The initial model thus obtained has 
the gravitational mass of M s = 1.72Mq and the radius of R s = 7.7km. This approximate 
neutron star is nearly the same as that used by Schinder et al (1988). 

In the dynamical simulation, we assume the same equation of state and use 100 
uniform baryon mass zones. In figures 8 and 9 the trajectories of all the mass zones are 
shown with the time variation of the Courant number. The time step is controlled so that 
all the dependent variables should not change more than 2% at each time step. Figure 8 
shows the mass element motions for the first 10msec. As can be seen, nothing occurs except 
for the very little oscillation during the first 1msec. The Courant number is shown with a 
curved line across the trajectories. It is found that the time step is increased monotonically. 
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This is because the matter velocities decrease monotonically after the first small transient. 
Figure 9 shows the mass trajectories for the same model but for a much longer time scale. 
The final time is 20sec, which should be compared with the dynamical time scale of the 
neutron star <^ 1msec. It is also noted this time corresponds to only 300 numerical time 
steps. The final Courant number exceeds 10 6 . These results imply that this dynamical 
calculation reproduced the hydrostatic configuration correctly and also guarantee that this 
code treats the general relativity appropriately. 



Finally shown are the results of the model calculations of the adiabatic core collapse. 
In these simulations, we assume the simplified equation of state which approximately 
reproduces the realistic simulations with detailed micro physics and neutrino transfer, and 
calculate the collapse of the iron core. The phenomenological equation of state adopted 
here is what Takahara & Sato (1982) used for their parametric research of the dynamics of 
core collapse. The total pressure p to t consists of the cold part p c and the thermal part p t as 
follows: 



where e t is the specific thermal energy density, and the adiabatic index F is assumed to be 
density dependent. Actually r = 4/3 for p b < 4 x 10 9 g/cm 3 and 10 12 < pb < 2.8 x 10 14 g/cm 3 , 
and T = 2.5 for p b > 2.8 x 10 14 g/cm 3 . The value of T for 4 x 10 9 < p b < 10 12 g/cm 3 is 
the model parameter. Similarly, the thermal stiffness 74 is varied from model to model. 
By adjusting these parameters, we can get both successful explosion models and stalled 
shock models. In fact, larger values help the shock propagation and tend to strengthen the 
explosion. The initial models are Woosley's 15Mq and 35Mq precollapse models ( |Woosley 
|1990|) . In the former model, the central 1.32Mq iron core is covered by 100 non-uniform 
mass zones in which the baryon mass of each cell is increased outward. T = 1.30 and 
7t = 1.30 are assumed. This model explodes in a prompt way. In the latter model, on the 
other hand, the central 2.0Mq is covered with 150 non-uniform meshes and T = 1.28 and 
7t = 1.25 are assumed. The shock is stalled inside the iron core for this model as shown 
below. 

Figures 10 and 11 show all the mass trajectories for the 15 Mq model with the time 
variation of the Courant number. Again I will show the results in two different time scales. 



2.7. The Adiabatic Collapse Calculation 



Ptot(Pb,£t) 
Pc(pb) 

Pt{Pb,£t) 



Pc(Pb) +Pt(pb,£t) 
Kpl 

(it ~ ±)pbSt , 



(61) 
(62) 
(63) 
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In figure 10, we can see the motions of mass elements until the shock reaches the iron core 
surface. As stated above, since the shock propagates the entire iron core without stagnation 
and some part of the outer core moves outward at the break out, this model is a successful 
prompt explosion model. It is also noted that the central part becomes almost hydrostatic. 
The data concerning the location of each mass element and the Courant number are saved 
every 10 time steps. As can be seen, the Courant number increases up to about 500 at the 
break out. Figure 11 shows the subsequent evolution until the calculation is terminated. 
The final time is 600msec and corresponds to 5000 numerical time steps. Since the mantle is 
not taken into consideration, the result after shock breaks out is not of physical significance. 
However, the calculation was continued in order to see how large the Courant number gets. 
As can be seen, the final Courant number exceeds 1000. One thing to be mentioned here 
is that the Courant number shown here is determined not at the shock wave but at the 
center of the proto-neutron star, which implies that we could overcome the severe Courant 
limit at the center without instability. This seems quite encouraging for the purpose of this 
project. It should be mentioned, however, that this is partly because there is no significant 
accretion on the proto-neutron star in the successful explosion model. Hence it is of greater 
importance to see the time evolution of the stalled shock model. 

In figure 12 we can see the motions of the representative mass elements for the 35Mq 
model until 50msec after the core bounce. In this case, the shock does not expel the 
accreting matter, that is, the shock is stalled. The inner core becomes nearly hydrostatic 
also in this case. However, since the matter is still falling onto the proto-neutron star 
continuously, the Courant number cannot become so large as the previous model as 
expected. The subsequent evolution is shown in figure 13. The final time is 900msec after 
the start of the simulation. The total number of numerical time steps is 8000. As can be 
seen, in this case the Courant number also becomes as large as 1000. Again this is a very 
encouraging result for this project. However, it should be mentioned that the resolution of 
the outer core is not sufficient to resolve the shock wave for the 35Mq model. This is not 
only because the variation of dependent variables is allowed up to 2% in this calculation 
(see the results of the relativistic shock tube problems), but also because matter continues 
to concentrate to the center. Hence we must implement some rezoning routine or we should 
use a larger number of zones. This is the problem to be solved in the next step. 

3. SUMMARY 

The implicit general relativistic hydrodynamical scheme is coded using the approximate 
Riemann solver and extensive numerical tests are accomplished. The scheme is very simple 
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since the formulae to evaluate the velocity and the pressure at the cell interface have such 
simple forms that we can easily calculate their derivatives with regard to the dependent 
variables. The numerical solutions obtained by this scheme show no oscillation after both 
shock waves and rarefaction waves unlike the schemes requiring an artificial viscosity, 
although a small erroneous jump appears after the contact surface. The calculations of the 
simple adiabatic collapse of massive stellar cores indicate this code is quite stable even if we 
take a much longer time step than the Courant time step, which fact is very encouraging 
for its application to the more realistic calculation of the late time explosion. However, 
it is also found that the time step should be adjusted so that the Courant number at the 
shock wave should be ~ 0.2. It is also probable some automatic rezoning routine should be 
implemented to resolve the shock wave sufficiently all through the calculation. 

The purpose of this project is to provide a tool to study the radiation hydrodynamics 
of supernova cores far beyond the dynamical time scale of the neutron star. We are now 
proceeding to the next step, in which the present hydrodynamic code is combined with a 
neutrino transfer code using the multi-energy group flux-limited diffusion scheme. However, 
the final goal of this project is to make a Boltzmann solver and combine it with this implicit 
hydrodynamical code thereby removing uncertainties regarding the neutrino transfer. 

In recent years, the importance of multi-dimensional effects on the supernova explosion 
has been pointed out and many detailed numerical simulations have been done. At present, 
it seems that the clues to the successful supernova explosion are such multi-dimensional 
effect and neutrinos. Although some effort has been made to do the simulation of multi- 
dimensional radiation transfer, it seems still far from satisfaction. Hence the one- dimensional 
simulations with detailed neutrino transfer and multi-dimensional simulations will play 
complementary roles to one another. It is also possible that one-dimensional simulation 
would serve as a guide to make some approximate treatment of the multi-dimensional 
neutrino transfer. 

I am grateful to S. Woosley for providing us with the details of his progenitor models. 
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Fig. 1. — The distribution of density, velocity, pressure and specific internal energy density at 
t = 0.675sec (400 time steps) are shown with the analytic solutions (solid lines) for 7 = 5/3. 

Fig. 2. — The distribution of density, velocity, pressure and specific internal energy density 
for the special relativistic shock tube problem described in the text at t — 3.03 x 10~ 10 sec 
(490 time steps) are shown with the analytic solutions (solid lines) for 7 = 5/3. 

Fig. 3. — The numerical solution for the same Riemann problem as Fig. 2 but with larger 
time step (the Courant number = 2.0). t = 2.95 x 10~ 10 sec (60 steps). 

Fig. 4. — The numerical solution for the same Riemann problem as Fig. 2 but with 
much larger time steps (the Courant number reaches ~ 8.0 at the end of the calculation). 
t = 2.80 x 10~ 10 sec (16 steps). 

Fig. 5. — The numerical solution for Sedov's point explosion problem described in the text. 
The time is 0.49sec and the step number is 10000. The solid lines are the exact solutions. 
7 = 5/3. 

Fig. 6. — The numerical solution for the self-similar collapse of a 7 = 4/3 polytrope. The 
distributions are shown every 100 steps. The times are 0.463, 0.700, 0.804, 0.881, 0.920, 
0.937 seconds. The solid lines are just connecting the data for the same time step. 

Fig. 7. — The trajectories of the representative mass elements for the Oppenheimer-Snyder 
dust collapse problem are shown with white circles. The solid lines are exact solutions. 
The baryon masses for the trajectories are 0.125, 0.375, 0.625, 0.875, 1.125, 1.375, 1.625, 
1.875M . 

Fig. 8. — The mass trajectories for the hydrostatic problem described in the text. The 
curved lines across the trajectories show the time variation of the Courant number. The 
evolution for only the first 10msec is shown. 

Fig. 9. — The same as Fig. 8 but for the whole calculation. 

Fig. 10. — The mass trajectories for the adiabatic collapse of the 15Mq model. The line 
across the trajectories show the time variation of the Courant number. The evolution only 
for the first 170msec is shown. 

Fig. 11. — The same as Fig. 10 but for the whole calculation. 

Fig. 12. — The mass trajectories for the adiabatic collapse of the 35Mq model. The line 
across the trajectories show the time variation of the Courant number. The evolution only 
for the first 400msec is shown. 
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Fig. 13. — The same as Fig. 12 but for the whole calculation. 



